function [ampphase_tab] = amp_phase(traces,... traces:irisFetch trace struct
    T,... %periods
    starttime,...
    endtime,...
    WaveletParameters,...
    ampphase_tab) %optional
% calculate amplitude and phase of specified periods at each channel
% Josh Crozier 2019

ampphaseplot = false;

if ~exist('ampphase_tab','var')
    ampphase_tab = table();
end

ExtendSignal = false;
VoicesPerOctave = 4;
padcycles = 3;

for tracein = 1:length(traces)
    for Tind = 1:length(T)
        periodrange = [T(Tind),T(Tind)*2];
    try
        column = strcat(traces(tracein).network,'_',...
            traces(tracein).station,'_',...
            traces(tracein).channel);
        
        local_starttime = max([starttime(Tind),traces(tracein).startTime]);
        startind = floor(seconds(local_starttime - traces(tracein).startTime)...
            *traces(tracein).sampleRate) + 1;
        local_endtime = min([endtime(Tind),traces(tracein).startTime...
            + seconds(traces(tracein).sampleCount/traces(tracein).sampleRate)]);
        endind = ceil(seconds(local_endtime - traces(tracein).startTime)...
            *traces(tracein).sampleRate) - 1;
        sampleCount = endind - startind + 1;
        
%         if sampleCount < 2
%            'debug pause' 
%         end
    
        %goertzel
        freq_indices = round((1./seconds(T(Tind)))/traces(tracein).sampleRate...
            *sampleCount) + 1;
        %taperwin = tukeywin(sampleCount,0.01);
        temp = goertzel(detrend(traces(tracein).data(startind:endind)),...
            freq_indices);
        temp = temp/(sampleCount/traces(tracein).sampleRate);
        ampphase_tab{Tind, strcat('amp_',column)} = abs(temp);
        ampphase_tab{Tind, strcat('phase_',column)} = angle(temp);
        
        if ampphaseplot
            %debug plot
            figure(5), clf
            subplot(4,1,1)
            tplot = (0:sampleCount-1)/traces(tracein).sampleRate;
            time = local_starttime + ...
                seconds(tplot);
            plot(time,detrend(traces(tracein).data(startind:endind)),...
                time,abs(temp)*cos(tplot*2*pi/seconds(T(Tind)) + angle(temp)))
            legend('full data','DFT component')
            title(strcat('T=',num2str(seconds(T(Tind)))))
            grid on
            xlim([min(time), max(time)])
        end
        
        %cwt to check for continuous phases
        %pad to minimize edge effects
        padlen = ceil(padcycles*seconds(T(Tind))*traces(tracein).sampleRate);
        padstartind = max([1, startind-padlen]);
        padendind = min([traces(tracein).sampleCount, endind+padlen]);
        [wt,~,~] = cwt(traces(tracein).data(padstartind:padendind),'morse',seconds(1/traces(tracein).sampleRate),...
            'ExtendSignal',ExtendSignal,...
            'periodrange',periodrange,...
            'VoicesPerOctave',VoicesPerOctave,...
            'WaveletParameters',WaveletParameters);
%         for jj = 1:size(wt,2)
%             wt(Twt>coi(jj),jj) = NaN;
%         end
        wta = angle(wt(1,1+padlen:end-padlen));
        if ampphaseplot
        	wtA = abs(wt(1,1+padlen:end-padlen));
        end
        unwrapwt = unwrap(wta);
        slope = 2*pi/seconds(T(Tind));
        t = (0:length(unwrapwt)-1)/traces(tracein).sampleRate;
        if ampphaseplot
            time = traces(tracein).startTime + ...
                seconds((padstartind + padlen)/traces(tracein).sampleRate) + ...
                seconds(t);
        end
        intercept = mean(unwrapwt)-slope*mean(t);
        fitline = intercept + t*slope;
        meanphaseerr = mean(abs(unwrapwt-fitline));
        ampphase_tab{Tind, strcat('errphase_',column)} = meanphaseerr;
        
        if ampphaseplot
            %debug plot
            subplot(4,1,2)
            plot(time,wtA,'k.')
            %title(strcat('T=',num2str(seconds(T(Tind)))))
            xlim([min(time), max(time)])
            ylabel('CWT Amplitude')
            grid on
            subplot(4,1,3)
            plot(time,wta,'k.',time,unwrapwt,'b.',time,fitline,'r')
            %title(strcat('T=',num2str(seconds(T(Tind)))))
            legend({'CWT Phase','Unwrapped CWT Phase','Phase for continuous signal'},...
                'Location','northwest')
            xlim([min(time), max(time)])
            ylabel('Angle (rad)')
            grid on
            subplot(4,1,4)
            plot(time,unwrapwt - fitline,'g')
            %title(strcat('T=',num2str(seconds(T(Tind)))))
            legend({'Deviation from continuous'},...
                'Location','northwest')
            xlim([min(time), max(time)])
            ylabel('Angle (rad)')
            grid on
            grid on
        end

    catch ME
        disp(getReport(ME))
        disp([startind, endind, length(traces(tracein).data)])
        disp(freq_indices)
    end
    end %end loop over T
end %end loop over stations

end

